#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import json
import base64
import math

import csv
import numpy as np
import matlab.engine
from PIL import Image

from models.data_model import Probe
from flask import make_response, abort
from loguru import logger

from webapi.matlab_engine import matlab_eng
from webapi.pai_image_analysis import PAIAnalysis
from webapi.pai_acquisition_data_manage import AcquisitionDataManage
import webapi.pai_global
# sys.path.append(os.getcwd())
# from utils import app_config

class MatlabMatEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.ndarray):
            return obj.tolist()
        return super().default(obj)

## 从采集的bin文件中提取信号
def extract_signal(bin_file, modality, selected_channels, valid_data_only, range):
    logger.debug("extract modality: {0} signal of channels: {1} from file: {2}. valid_data_only :{3} , range: {4}.", modality, selected_channels, bin_file, valid_data_only, range)
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    pai_analysis.setCurrentAcquisitionFile(bin_file)
    currentAcqPath = pai_analysis.getCurrentAcquisitionPath()
    if currentAcqPath is None or not currentAcqPath.strip() :
        abort("CurrentAcquisitionPath is None! Call GET /acquisition/{acquisition_id} first to set CurrentAcquisitionPath.", 500)
    
    bin_file_path = currentAcqPath + "/" + bin_file
    acq_parameter_file_path = currentAcqPath + "/" + "acquisition_parameter.csv"

    matlab_scripts_spectrum = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_spectrum)

    modality_mode = {"US": 0, "PA": 1}.get(modality, -1)
    valid_flag = 1 if valid_data_only else 0

    py_channels = []
    for channel in selected_channels:
        str_num = float(channel[2:].strip())
        py_channels.append(str_num)
    matlab_channels = matlab.double(py_channels)
    matlab_range = matlab.double(range)

    logger.debug("call matlab function:PAUSSignalExtract. bin_file_path:{0}, modality_mode:{1}, matlab_channels:{2}, valid_flag:{3},matlab_range:{4}", bin_file_path, modality_mode, matlab_channels, valid_flag, matlab_range)
    extracted_signal = matlab_eng.PAUSSignalExtract(bin_file_path, 
                                                    acq_parameter_file_path,
                                                    float(modality_mode), 
                                                    matlab_channels, 
                                                    float(valid_flag), 
                                                    matlab_range)
    signal = np.array(extracted_signal["signals"]._data)
    signal = signal.transpose()
    extracted_signal["signals"] = signal

    XAxisTickLabel = np.array(extracted_signal["XAxisTickLabel"]._data)
    extracted_signal["XAxisTickLabel"] = XAxisTickLabel

    YAxisTickLabel = np.array(extracted_signal["YAxisTickLabel"]._data)
    extracted_signal["YAxisTickLabel"] = YAxisTickLabel

    selectedchannels = np.array(extracted_signal["selectedchannels"]._data)
    extracted_signal["selectedchannels"] = selectedchannels

    return json.dumps(extracted_signal, cls = MatlabMatEncoder), 200

## TBD 后续可优化，不用每次都读文件取数据
def signal_spectrum(bin_file, modality, selected_channels, valid_data_only, range, method):
    logger.debug("signal madality: {0} signal of channels: {1} from file: {2}. valid_data_only :{3} , range: {4}, method: {5}.", modality, selected_channels, bin_file, valid_data_only, range, method)
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    pai_analysis.setCurrentAcquisitionFile(bin_file)
    currentAcqPath = pai_analysis.getCurrentAcquisitionPath()
    if currentAcqPath is None or not currentAcqPath.strip() :
        abort("CurrentAcquisitionPath is None! Call GET /acquisition/{acquisition_id} first to set CurrentAcquisitionPath.", 500)
    
    bin_file_path = currentAcqPath + "/" + bin_file
    acq_parameter_file_path = currentAcqPath + "/" + "acquisition_parameter.csv"

    matlab_scripts_spectrum = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_spectrum)

    modality_mode = {"US": 0, "PA": 1}.get(modality, -1)
    valid_flag = 1 if valid_data_only else 0

    py_channels = []
    for channel in selected_channels:
        str_num = float(channel[2:].strip())
        py_channels.append(str_num)
    matlab_channels = matlab.double(py_channels)
    matlab_range = matlab.double(range)

    logger.debug("call matlab function:PAUSSignalSpectrumFromBin. bin_file_path:{0}, modality_mode:{1}, matlab_channels:{2}, valid_flag:{3},matlab_range:{4},method:{5}.", bin_file_path, modality_mode, matlab_channels, valid_flag, matlab_range, method)
    extracted_signal = matlab_eng.PAUSSignalSpectrumFromBin(bin_file_path, 
                                                            acq_parameter_file_path,
                                                            float(modality_mode), 
                                                            matlab_channels, 
                                                            float(valid_flag), 
                                                            matlab_range,
                                                            method)
    signal = np.array(extracted_signal["signals"]._data)
    signal = signal.transpose()
    extracted_signal["signals"] = signal

    XAxisTickLabel = np.array(extracted_signal["XAxisTickLabel"]._data)
    extracted_signal["XAxisTickLabel"] = XAxisTickLabel
    #print(f"XAxisTickLabel:", extracted_signal["XAxisTickLabel"])
    #extracted_signal["XAxisLabel"]

    YAxisTickLabel = np.array(extracted_signal["YAxisTickLabel"]._data)
    extracted_signal["YAxisTickLabel"] = YAxisTickLabel
    #print(f"YAxisTickLabel:", extracted_signal["YAxisTickLabel"])

    selectedchannels = np.array(extracted_signal["selectedchannels"]._data)
    extracted_signal["selectedchannels"] = selectedchannels
    #print(extracted_signal["selectedchannels"])

    return json.dumps(extracted_signal, cls = MatlabMatEncoder), 200

## get supported reconstruction methods
def get_supported_recon_methods():
    logger.debug("start get_supported_recon_methods.")
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    recon_methods_list = pai_analysis.getsupportedReconMethods()
    return recon_methods_list, 200

## 图像重建
def image_reconstruction(bin_file, modality, method, log_compress_parameter, sound_speed, generate_crop_roi):
    logger.debug("start image_reconstruction. bin file: {0}, modality: {1}, reconstruction medthod {2}, log_compress_parameter :{3} , sound_speed: {4}, generate_crop_roi: {5}", bin_file, modality, method, log_compress_parameter, sound_speed, generate_crop_roi)
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    pai_analysis.setCurrentAcquisitionFile(bin_file)
    currentAcqPath = pai_analysis.getCurrentAcquisitionPath()
    if currentAcqPath is None or not currentAcqPath.strip() :
        abort("CurrentAcquisitionPath is None! Call GET /acquisition/{acquisition_id} first to set CurrentAcquisitionPath.", 500)
    
    bin_file_path = currentAcqPath + "/" + bin_file
    acq_parameter_file_path = currentAcqPath + "/" + "acquisition_parameter.csv"

    # 获取探头的参数
    with open(acq_parameter_file_path, "r", encoding="utf-8") as f:
        csv_reader = csv.DictReader(f)
        for row in csv_reader:
            probe_name = row['ProbeName']

    acq_data_manager = webapi.pai_global.get_value("AcquisitionDataManage")
    probe = acq_data_manager.getProbeInfo(probe_name)

    probe_params = {"Nelements":float(probe.Nelements), 
                    "fc":float(probe.fc), 
                    "pitch":float(probe.pitch), 
                    "kerf":float(probe.kerf), 
                    "width":float(probe.width), 
                    "bandwidth":float(probe.bandwidth), 
                    "radius":float(probe.radius), 
                    "focus":float(probe.focus), 
                    "height":float(probe.height),
                    "name":probe.name}
    #调用matlab函数
    matlab_scripts_path = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_path)
    matlab_must_path = matlab_scripts_path + "/MUST"
    matlab_eng.addpath(matlab_must_path)

    modality_mode = {"US": 0, "PA": 1}.get(modality, -1)
    
    save_path = currentAcqPath
    logger.debug("call matlab function: USPArecon. bin_file_path:{0}, probe_params:{1}, modality_mode:{2}, reconstruction method:{3} , log_compress_parameter:{4},  sound_speed{5}:", bin_file_path, probe_params, modality_mode, method, log_compress_parameter, sound_speed)
    recon_image = matlab_eng.USPArecon(bin_file_path, 
                                       acq_parameter_file_path,
                                       save_path,
                                       probe_params,
                                       modality_mode,
                                       method,
                                       float(log_compress_parameter),
                                       float(sound_speed))
    #适配输出结果到前端
    image_data = np.array(recon_image["image"]._data)
    dimension = np.array(recon_image["dimension"]._data)
    recon_image["dimension"] = dimension

    pixel_spacing = np.array(recon_image["pixel_spacing"]._data)
    recon_image["pixel_spacing"] = pixel_spacing

    Colormap = np.array(recon_image["Colormap"]._data) * 255
    Colormap = Colormap.astype(np.uint8)
    color = Colormap.reshape(256, 3)
    palette = np.zeros((256, 4), dtype=np.uint8)
    palette[:, 0] = color[:, 2]
    palette[:, 1] = color[:, 1]
    palette[:, 2] = color[:, 0]
    palette[:, 3] = 0
    recon_image["Colormap"] = Colormap

    XTick = np.array(recon_image["XTick"]._data)
    recon_image["XTick"] = XTick

    YTick = np.array(recon_image["YTick"]._data)
    recon_image["YTick"] = YTick

    bmp_data = pai_analysis.generateBMPWithPaletteImage(image_data, dimension, palette)
    ## code for testing
    # with open('pixeldata.dat', 'wb') as f:
    #     f.write(bmp_data)
    #pai_analysis.showImage(bmp_data)
    curAcqID = pai_analysis.getCurrentAcquisitionID()
    pai_analysis.setImage(curAcqID, modality, bmp_data, {"pixel_spacing":pixel_spacing, "XTick":XTick, "YTick":YTick})
    
    recon_image["image"] = str(base64.b64encode(bmp_data))
    ## code for testing
    # with open('output2.dat', 'wb') as f:
    #     f.write(bmp_data)

    if generate_crop_roi:
        crop_box = pai_analysis.getCropbox()
        crop_image = pai_analysis.cropImage(bmp_data, crop_box)
        crop_bmp_image = pai_analysis.image2BMPBytes(crop_image)
        xsmall = min(crop_box[0], crop_box[2])
        xbig = max(crop_box[0], crop_box[2])
        xtick1 = math.floor(xsmall * pixel_spacing[0])
        xtick2 = math.ceil(xbig * pixel_spacing[0])
        xdelta = round((xtick2 - xtick1) / 4)
        XTick = np.array([xtick1, xtick1+xdelta, xtick1+2*xdelta, xtick1+3*xdelta, xtick1+4*xdelta])

        ysmall = min(crop_box[1], crop_box[3])
        ybig = max(crop_box[1], crop_box[3])
        ytick1 = math.floor(ysmall * pixel_spacing[1])
        ytick2 = math.ceil(ybig * pixel_spacing[1])
        ydelta = round((ytick2 - ytick1) / 4)
        YTick = np.array([ytick1, ytick1+ydelta, ytick1+2*ydelta, ytick1+3*ydelta, ytick1+4*ydelta])

        pai_analysis.setCropImage(curAcqID, modality, crop_bmp_image, {"pixel_spacing":pixel_spacing, "XTick":XTick, "YTick":YTick})
        recon_image["crop_roi_image"] = str(base64.b64encode(crop_bmp_image))
        recon_image["crop_image_YTick"] = YTick
        recon_image["crop_image_XTick"] = XTick
        #recon_image["crop_image_XTick"] = np.array([25, 30, 35])

    logger.debug("end image_reconstruction.")
    return json.dumps(recon_image, cls = MatlabMatEncoder), 200

## 图像融合
def paus_image_fusion(bin_file, generate_crop_roi):
    logger.debug("start paus_image_fusion input parameter: bin_file - {0}, generate_crop_roi - {1}", bin_file, generate_crop_roi)
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    current_bin_file = pai_analysis.getCurrentAcquisitionFile()
    if current_bin_file != bin_file:
        return "the input bin file is not consistent with the current bin file stored in the server.", 400

    curAcqID = pai_analysis.getCurrentAcquisitionID()
    fusion_image = pai_analysis.pausFusion(curAcqID, 0.5)
    if "XTick" not in fusion_image.info or "YTick" not in fusion_image.info:
        abort("lack of XTick or YTick in fusion_image.info. please add XTick or YTick info in image_reconstruction.")

    fusion_bmp_image = pai_analysis.image2BMPBytes(fusion_image)
    fusion_image = {"image": str(base64.b64encode(fusion_bmp_image)),
                    "XTick": fusion_image.info["XTick"],
                    "YTick": fusion_image.info["YTick"]}

    if generate_crop_roi:
        crop_box = pai_analysis.getCropbox()
        crop_fusion_image = pai_analysis.cropImage(fusion_bmp_image, crop_box)

        modality = "US"
        acq_id_list = [curAcqID]
        crop_us_image_list = pai_analysis.getCropImageList(modality, acq_id_list)
        crop_us_image = crop_us_image_list[0]
        if "XTick" not in crop_us_image.info or "YTick" not in crop_us_image.info:
            abort("lack of XTick or YTick in crop_us_image.info. please add XTick or YTick info in image_reconstruction.")

        crop_fusion_bmp_image = pai_analysis.image2BMPBytes(crop_fusion_image)
        fusion_image.update({"crop_roi_image": str(base64.b64encode(crop_fusion_bmp_image)),
                             "crop_image_YTick": crop_us_image.info["YTick"],
                             "crop_image_XTick": crop_us_image.info["XTick"]})

    # fusion_image["dimension"]
    # fusion_image["pixel_spacing"]
    # fusion_image["XTick"]
    # fusion_image["YTick"]

    logger.debug("end paus_image_fusion.")
    return json.dumps(fusion_image, cls = MatlabMatEncoder), 200

## 复合平面波仿真
def simulation_compound_plane_wave():
    logger.debug("start simulation_compound_plane_wave。")
    matlab_scripts_path = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_path)

    matlab_must_path = matlab_scripts_path + "/MUST"
    matlab_eng.addpath(matlab_must_path)

    logger.debug("call matlab function:CompoundPlaneWaveSimulation. ")
    [angle0_result, angle6_result, compound_result, comparison_result] = matlab_eng.CompoundPlaneWaveSimulation(nargout = 4)
    
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")

    angle0_image_data = np.array(angle0_result["image"]._data)
    dimension = np.array(angle0_result["dimension"]._data)
    angle0_result["dimension"] = dimension
    pixel_spacing = np.array(angle0_result["pixel_spacing"]._data)
    angle0_result["pixel_spacing"] = pixel_spacing
    xposition = np.array(angle0_result["Xposition"]._data)
    angle0_result["Xposition"] = xposition
    zposition = np.array(angle0_result["Zposition"]._data)
    angle0_result["Zposition"] = zposition

    palette = np.zeros((256, 4), dtype=np.uint8)
    palette[:, 0] = np.arange(256, dtype=np.uint8)
    palette[:, 1] = np.arange(256, dtype=np.uint8)
    palette[:, 2] = np.arange(256, dtype=np.uint8)
    palette[:, 3] = 0
    angle0_result["Colormap"] = palette

    bmp_data = pai_analysis.generateBMPWithPaletteImage(angle0_image_data, dimension, palette)

    #bmp_data = pai_analysis.generateBMPGrayImage(angle1_image_data, dimension)
    angle0_result["image"] = str(base64.b64encode(bmp_data))
    # # code for testing
    # with open('output2.dat', 'wb') as f:
    #     f.write(bmp_data)

    angle6_image_data = np.array(angle6_result["image"]._data)
    dimension = np.array(angle6_result["dimension"]._data)
    angle6_result["dimension"] = dimension
    pixel_spacing = np.array(angle6_result["pixel_spacing"]._data)
    angle6_result["pixel_spacing"] = pixel_spacing
    xposition = np.array(angle6_result["Xposition"]._data)
    angle6_result["Xposition"] = xposition
    zposition = np.array(angle6_result["Zposition"]._data)
    angle6_result["Zposition"] = zposition
    bmp_data = pai_analysis.generateBMPWithPaletteImage(angle6_image_data, dimension, palette)
    angle6_result["image"] = str(base64.b64encode(bmp_data))

    compound_image_data = np.array(compound_result["image"]._data)
    dimension = np.array(compound_result["dimension"]._data)
    compound_result["dimension"] = dimension
    pixel_spacing = np.array(compound_result["pixel_spacing"]._data)
    compound_result["pixel_spacing"] = pixel_spacing
    xposition = np.array(compound_result["Xposition"]._data)
    compound_result["Xposition"] = xposition
    zposition = np.array(compound_result["Zposition"]._data)
    compound_result["Zposition"] = zposition
    bmp_data = pai_analysis.generateBMPWithPaletteImage(compound_image_data, dimension, palette)
    compound_result["image"] = str(base64.b64encode(bmp_data))

    signal = np.array(comparison_result["OneAngle"]._data)
    signal = signal.transpose()
    comparison_result["OneAngle"] = signal
    signal = np.array(comparison_result["compound"]._data)
    signal = signal.transpose()
    comparison_result["compound"] = signal
    xAxisTickLabel = np.array(comparison_result["XAxisTickLabel"]._data)
    comparison_result["XAxisTickLabel"] = xAxisTickLabel
    yAxisTickLabel = np.array(comparison_result["YAxisTickLabel"]._data)
    comparison_result["YAxisTickLabel"] = yAxisTickLabel

    compound_plane_wave_simulation = {"angle0_result": json.dumps(angle0_result, cls = MatlabMatEncoder),
                       "angle6_result": json.dumps(angle6_result, cls = MatlabMatEncoder),
                       "compound_result": json.dumps(compound_result, cls = MatlabMatEncoder),
                       "comparison_result": json.dumps(comparison_result, cls = MatlabMatEncoder)}
    
    logger.debug("end simulation_compound_plane_wave。.")
    return json.dumps(compound_plane_wave_simulation, cls = MatlabMatEncoder), 200

## 全波反演仿真
def simulation_full_waveform_inversion(anatomic_part_name, inversion_parameter):
    logger.debug("start simulation_full_waveform_inversion input anatomic_part_name: {0}, inversion_parameter:{1}", anatomic_part_name, inversion_parameter)
    
    if anatomic_part_name not in ["tibia_and_fibula", "tibia"]:
        return "input parameter anatomic_part_name is invalid. ", 501
    if inversion_parameter not in ["density", "velocity"]:
        return "input parameter anatomic_part_name is invalid. ", 501

    matlab_scripts_path = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_path)

    logger.debug("call matlab function:FullWaveformInversionSimulation. anatomic_part_name:{0}, inversion_parameter:{1}",anatomic_part_name, inversion_parameter)
    bmp_data = matlab_eng.FullWaveformInversionSimulation(anatomic_part_name, inversion_parameter)
    # ## code for testing
    # with open('output2.dat', 'wb') as f:
    #     f.write(bytes(bmp_data._data))
    simulation_result = {"image": str(base64.b64encode(bytes(bmp_data._data)))}
    logger.debug("end simulation_full_waveform_inversion.")
    return simulation_result, 200

## 分子成像，解混不同化合物的成像。 TBD 暂时为mock函数，后续需根据实际情况重写。
def pa_molecular_imaging_unmixing(acquisition_id_list, for_fe_test):
    logger.debug("start molecular_imaging_unmixing. input acquisition_id_list: {0}. for_fe_test:{1}", acquisition_id_list, for_fe_test)
    
    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")
    if for_fe_test:
        with open("./matlab_scripts/const_file/mineral_result.bmp", "rb") as fid:
            mineral_result_bmp = fid.read()
        with open("./matlab_scripts/const_file/lipid_result.bmp", "rb") as fid:
            lipid_result_bmp = fid.read()
        with open("./matlab_scripts/const_file/collagen_result.bmp", "rb") as fid:
            collagen_result_bmp = fid.read()
        with open("./matlab_scripts/const_file/Hb_result.bmp", "rb") as fid:
            Hb_result_bmp = fid.read()
        with open("./matlab_scripts/const_file/HHb_result.bmp", "rb") as fid:
            HHb_result_bmp = fid.read()
        with open("./matlab_scripts/const_file/BOS_result.bmp", "rb") as fid:
            BOS_result_bmp = fid.read()

        unmixing_result = {"mineral": str(base64.b64encode(mineral_result_bmp)),
                       "lipid": str(base64.b64encode(lipid_result_bmp)),
                       "collagen": str(base64.b64encode(collagen_result_bmp)),
                       "Hb": str(base64.b64encode(Hb_result_bmp)),
                       "HHb": str(base64.b64encode(HHb_result_bmp)),
                       "BOS": str(base64.b64encode(BOS_result_bmp))}
        return unmixing_result, 200

    if len(acquisition_id_list) < 7:
        return "at least 6 images of different laser wavelengths are needed. ", 500

    matlab_scripts_path = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_path)

    pa_image_list = pai_analysis.getCropImageList("PA", acquisition_id_list)
    us_acq_id_list = [acquisition_id_list[0], acquisition_id_list[-1]]
    us_image_list = pai_analysis.getCropImageList("US", us_acq_id_list)
    us_image = us_image_list[0]

    pa_image_pixel_list = []
    for img in pa_image_list:
        img_pixel_bytes = pai_analysis.extractImagePixelBytes(img)
        pa_image_pixel_list.append(matlab.uint8(img_pixel_bytes))
    
    us_pixel_bytes = matlab.uint8(pai_analysis.extractImagePixelBytes(us_image))
    image_width = us_image.width
    image_height = us_image.height

    wavelength_list = []
    acq_data_manage = webapi.pai_global.get_value("AcquisitionDataManage")
    for acq_id in acquisition_id_list:
        acquisition = acq_data_manage.getAcquisitionInfoByAcquisitionID(acq_id)
        wavelength_list.append(acquisition.laser_wavelength)
    wavelength_list = matlab.double(wavelength_list)

    logger.debug("call matlab function: PAMolecularImagingUnmixing. count of PA images:{0}, image width:{1}, image height:{2}, input laser wavelength list:{3}", len(pa_image_pixel_list), image_width, image_height, wavelength_list)
    [mineral_result, lipid_result, collagen_result, Hb_result, HHb_result, BOS_result] = matlab_eng.PAMolecularImagingUnmixing(pa_image_pixel_list, us_pixel_bytes, float(image_width), float(image_height), wavelength_list, nargout = 6)
    mineral_result_bmp = pai_analysis.generateBMPImage(np.array(mineral_result["image"]._data), np.array(mineral_result["dimension"]._data))
    # ## code for testing
    # with open('output.dat', 'wb') as f:
    #     f.write(mineral_result_bmp)
    lipid_result_bmp = pai_analysis.generateBMPImage(np.array(lipid_result["image"]._data), np.array(lipid_result["dimension"]._data))
    collagen_result_bmp = pai_analysis.generateBMPImage(np.array(collagen_result["image"]._data), np.array(collagen_result["dimension"]._data))
    Hb_result_bmp = pai_analysis.generateBMPImage(np.array(Hb_result["image"]._data), np.array(Hb_result["dimension"]._data))
    HHb_result_bmp = pai_analysis.generateBMPImage(np.array(HHb_result["image"]._data), np.array(HHb_result["dimension"]._data))
    BOS_result_bmp = pai_analysis.generateBMPImage(np.array(BOS_result["image"]._data), np.array(BOS_result["dimension"]._data))

    unmixing_result = {"mineral": str(base64.b64encode(mineral_result_bmp)),
                       "lipid": str(base64.b64encode(lipid_result_bmp)),
                       "collagen": str(base64.b64encode(collagen_result_bmp)),
                       "Hb": str(base64.b64encode(Hb_result_bmp)),
                       "HHb": str(base64.b64encode(HHb_result_bmp)),
                       "BOS": str(base64.b64encode(BOS_result_bmp))}
    logger.debug("end molecular_imaging_unmixing.")
    return unmixing_result, 200

## 基于重建的光声或超声图像测量骨小梁感兴趣区的平均厚度
def trabecular_thickness_measurement(acquisition_id, modality, use_full_image, roi_type, point_list):
    logger.debug("start trabecular_thickness_measurement. acquisition_id:{0}, modality:{1}, use_full_image:{2}, roi_type:{3}, point_list{4}", acquisition_id, modality, use_full_image, roi_type, point_list)
    
    if roi_type != "rectangle":
        return "only roi_type rectangle is supported now. ", 501

    matlab_scripts_path = os.getcwd() + "/matlab_scripts/image_signal_analysis"
    matlab_eng.cd(matlab_scripts_path)

    pai_analysis = webapi.pai_global.get_value("PAIAnalysis")

    acq_id_list = [acquisition_id]
    
    if use_full_image:
        image_list = pai_analysis.getImageList(modality, acq_id_list)
    else:
        image_list = pai_analysis.getCropImageList(modality, acq_id_list)
    
    if len(image_list) == 0:
        return "There is no image matched with the input acquisition_id", 501
    
    image = image_list[0]
    logger.debug("trabecular_thickness_measurement image info:{0}", image.info)
    #print(image.info)
    if image.info["pixel_spacing"] != None:
        pixel_spacing = float(image.info["pixel_spacing"][0])
    else:
        return "there is no pixel_spacing in PIL image.info. Cannot calculate trabecular_thickness.", 501

    roi_image = pai_analysis.extractROIImage(image, roi_type, point_list)
    if not roi_image:
        return "cannot extract roi sub image according to the imput point list.", 501

    pixel_bytes = matlab.uint8(pai_analysis.extractImagePixelBytes(roi_image))
    image_width = roi_image.width
    image_height = roi_image.height

    logger.debug("call matlab function: TrabecularThicknessMeasurement. image width:{0}, image height:{1}, pixel_spacing:{2}", image_width, image_height, pixel_spacing)
    thickness = matlab_eng.TrabecularThicknessMeasurement(pixel_bytes, float(image_width), float(image_height), float(pixel_spacing))
    trabecular = {"trabecular thickhness":thickness, "unit": "mm"}
    return trabecular, 200


if __name__ == "__main__":
    pass